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AN  IMPROVED  PROCEDURE  FOR  ADAPTIVE  ESTIMATION 
OF  MANEUVERING  IN  A  SHIP  TRACKING  ALGORITHM 


1 .  INTRODUCTION 

A  single-target  ship  tracking  algorithm  has  been  developed  in  Refs.  [1]  - 
C33,  and  is  also  incorporated  as  a  component  of  the  multi-target  tracker 
described  in  Ref.  [4],  This  algorithm  operates  on  ship  location  reports  which 
may  occur  sporadically  in  time,  and  adaptively  estimates  the  intensity  of  the 
ship's  maneuvering  behavior  from  the  report  data  rather  than  accepting  it  as  a 
user-specified  parameter.  Maneuvering  here  is  defined  as  deviations  from  a 
steady  course,  which  is  also  estimated  from  the  data.  Ref.  [33  also  describes 
a  variant  of  this  algorithm  which  does  not  distinguish  between  down-track  and 
cross-track  maneuvering,  but  which  can  make  more  effective  use  of  bearing  - 
only  reports,  as  opposed  to  reports  of  ship  position  which  are  meaningfully 
localized  in  two  dimensions. 

Another  procedure  has  subsequently  been  devised  for  this  maneuvering  esti¬ 
mation  which  significantly  improves  the  performance  of  both  tracking  algorithm 
variants.  This  new  procedure  is  the  subject  of  this  report,  and  differs 
basically  by  the  use  of  weighted  averages  of  certain  statistics,  rather  than 
simple  averages,  to  form  estimates  of  maneuvering  intensity.  The  reliability 
of  these  statistics  as  measurements  of  maneuvering  intensity  varies  with  the 
spacing  of  adjacent  reports.  The  new  procedure  gives  relatively  lower  weight 
to  the  less  reliable  statistics  in  the  case  of  unevenly  spaced  reports, 
resulting  in  more  reliable  maneuvering  intensity  estimates.  Hence,  the 
resulting  improvement  in  tracking  performance  is  most  pronounced  when  the 
report  occurrances  are  highly  irregular.  However,  there  seems  to  be  3ome 
improvement  even  in  the  case  of  regularly  spaced  reports,  as  explained  in 
Section  6  of  this  report. 

2.  PURPOSE  AND  ROLE  OF  IMPROVEMENT 

For  planar  motion,  the  ship  tracking  algorithm  of  Refs.  [1]  -  [33  is  based 
on  the  idea  of  approximating  the  motion  of  the  (single)  ship  being  tracked  as 
the  vector  sum  of  a  constant  (average)  velocity  and  a  two-dimensional  (random) 
Brownian  motion.  The  Brownian  motion's  Intensity  parameter,  a  2x2  matrix  Q  in 
rectangular  coordinates,  is  selected  to  correspond  to  the  extent  of  maneuver¬ 
ing  performed  by  the  ship  with  respect  to  a  constant-speed,  constant-heading 
Manuscript  submitted  March  2, 1982. 
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course  (the  average  velocity).  Both  the  average  velocity  and  the  maneuvering 
intensity  matrix  Q  are  treated  as  constant  parameters  to  be  estimated  from  the 
observed  input  data.  This  input  consists  of  a  time-ordered  series  of  reports 
of  ship  location,  not  necessarily  evenly  spaced  in  time,  such  that  each  report 
specifies  a  time,  the  observed  ship  position  at  that  time,  and  the  (2x2) 
covariance  matrix  (or,  equivalently,  a  containment  ellipse)  for  the  error  in 
this  observed  position.  An  exception,  however,  is  a  modified  tracking  proced¬ 
ure  described  in  Ref.  [3],  which  allows  some  (or  all)  of  these  reports  to 
specify  only  the  observed  ship  bearing  from  some  given  point  at  that  time,  and 
the  variance  of  the  bearing  observation  error. 

In  either  case,  the  tracking  algorithm  operates  recursively  in  time, 
basically  by  propagating  the  ship  track  forward  between  observations  by  dead 
reckoning  and  updating  it  whenever  a  new  report  is  received.  The  ship 
position  and  average  velocity  are  treated  as  a  four-component  "state  vector," 
for  which  a  current  conditional  mean  and  covariance  matrix  are  generated  by  a 
standard  Kalman  filter.  Another  recursive  procedure  is  used  to  estimate  the 
maneuvering  intensity  parameter  Q  from  the  "innovations"  of  the  Kalman  filter. 
These  estimates  are  then  used  as  "driving  noise"  parameters  in  the  Kalman 
filter  adaptively  to  modify  its  subsequent  operation.  The  details  of  this 
process  are  described  in  Refs.  [1]  -  [3],  but  can  be  summarized  as  the 
following  recurring  sequence  of  basic  steps: 

1.  Upon  the  receipt  of  a  new  report  of  observed  ship  location,  propagating 
the  conditional  probability  distribution  of  the  ship's  "state  vector" 
(position  and  velocity),  given  all  previous  data,  to  the  current  time. 

Thi3  distribution  (conditioned  on  the  same  data)  is  generally  available 
for  the  3tate  at  some  earlier  time,  and  is  propagated  from  that  time  as  if 
the  value  of  Q  estimated  then  were  a  precisely  known  parameter.  The  state 
distribution  is  treated  as  if  it  were  (4-variate)  Normal,  so  this  step 
just  amounts  to  propagating  its  mean  and  covariance  matrix  (14  independent 
components)  to  the  current  time  with  the  standard  Kalman  filter  for  this 
value  of  Q  and  the  motion  model  described  above. 

2.  Updating  the  state  probability  distribution  with  the  new  observation, 
again  using  the  standard  Kalman  filter  for  this  case.  This  gives  the 
conditional  state  distribution  (specified  by  the  mean  and  covariance 
matrix)  given  the  current  observation  as  well. 
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3.  Updating  the  estimate  of  Q  using  the  innovations  (observed  minus 
propagated  position)  from  step  2. 

4.  Adjustment  of  some  of  the  parameters  of  the  updated  position-velocity 
distribution  from  step  2  to  compensate  for  the  fact  that  the  value  of  Q 
used  in  step  1  was  (in  general)  different  from  the  one  just  estimated  in 
step  3. 

Steps  1  and  4  are  not  performed  for  the  first  report.  In  this  case,  the  track 
is  initiated  with  an  estimated  Q  of  zero  and  a  user-specified  zero-mean 
circular  Normal  distribution  for  the  average  velocity.  Also,  step  1  can  be 
performed  to  project  a  ship  location  distribution  to  a  time  at  which  there  is 
no  observation.  This  just  amounts  to  an  extended  form  of  dead  reckoning,  in 
which  an  entire  containment  ellipse  (representing  the  Normal  distribution  of 
the  position  components  of  the  state  vector)  is  propagated,  not  just  a  most 
likely  position  (the  center  of  this  ellipse).  In  step  3.  it  is  assumed  that 
the  Brownian  maneuvering  motion  consists  of  statistically  independent  compon¬ 
ents  parallel  and  perpendicular  to  the  average  velocity,  so  only  a  down-track 
and  a  cross-track  maneuvering  intensity  are  estimated,  and  the  resulting  Q 
matrix  is  always  diagonal  when  tranformed  for  coordinates  aligned  with  the 
currently  estimated  average  velocity  vector.  The  intensiities  of  these  two 
Brownian  motion  components  are  further  restricted  to  be  identical  (isotropic 
maneuvering)  in  the  variant  of  Ref.  [33,  which  allows  the  inclusion  of 
bearing-only  data. 

The  improvement  reported  here  is  in  the  procedure  for  estimating  the 
maneuvering  intensity  parameter  Q  in  step  3.  The  advantage  of  thi3  new  pro¬ 
cedure  occurs  chiefly  when  the  observation  times  are  quite  unevenly  spaced. 

The  difficulty  with  the  previous  method  is  that,  in  both  variants,  equal 
weight  is  given  to  each  member  of  the  Kalman  filter’s  innovation  sequence 
(from  step  2)  in  constructing  the  estimate  of  Q.  However,  the  precision  of 
these  innovation  statistics  as  measure  of  the  Q  components  varies  approximate¬ 
ly  as  the  inverse  of  the  elapsed  time  since  the  last  report,  so  their  use  in 
this  way  for  unusually  closely  spaced  reports  often  results  in  a  less  reliable 
estimate  of  Q  than  ignoring  them  entirely.  The  new  estimation  procedure  is 
devised  so  that  these  innovation  statistics  are  weighted  according  to  their 
reliability  as  measures  of  the  Q-components  being  estimated,  thereby  avoiding 
this  problem.  As  with  the  former  procedure,  variants  can  be  developed  for 
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both  the  case  of  independent  down-track  and  cross-track  maneuvering  components 
and  the  case  of  isotropic  maneuvering  with  bearing-only  reports,  as  described 
in  detail  below. 

Another  possible  benefit  of  the  new  Q-estimation  procedure  described  here 
is  that  it  seems  to  cause  the  overall  tracking  algorithm  to  create  somewhat 
larger  containment  ellipses  for  projected  ship  locations  than  the  former  pro¬ 
cedure,  even  for  evenly  spaced  observation  times  (see  Figs.  2  and  3).  This 
would  correct  a  reported  tendency  of  the  former  method  to  result  in  erroneous 
ly  small  containment  ellipses  when  used  on  realistic  ship  tracking  data. 


Fig.  1  —  Outline  of  single-target  ship  tracking  algorithm 
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3.  RATIONALE  FOR  ONE-DIMENSIONAL  SHIP  MOTION 

The  new  procedure  for  estimating  the  ship's  maneuvering  intensity 
parameter  Q  is  based  on  a  consideration  of  the  corresponding  one-dimensional 
case,  in  which  the  ship's  motion  is  described  by  the  (scalar)  equation 

X  s  u  +  w 


and  the  position  observations  at  discrete  time  t^  are 

zi  8  x<ti>  "*■  ni*  1  =  1 . 

where  the  n^  are  independent  zero-mean  Normal  random  variables  with  known 

variances  r^,  u  is  a  constant  but  a-priori  unknown  average  velocity,  and  w  is 

a  Normal  white  noise  process  with  constant,  but  a-priori  unknown  variance 

parameter  q.  Tracking  begins  immediately  after  the  initial  observation  zq  at 

time  t  ,  at  which  point  the  two-vector 
o 


is  regarded  as  having  a  bivariate  Normal  probability  distribution  with  mean 


and  covariance  matrix 


EH'S . 


2 

where  x  is  some  a-priori  specified  ship  speed  variance. 

Now  consider  the  approximation  that,  at  observation  time  t  ,  the 
conditional  probability  densities  of  q  and  of  the  vector  jx(t^,  given 


the  observations  z  z, ,  are  statistically  independent  such  that, 

o  1 


q  —  Normal  (q,  s) 
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and 


Normal 


For  the  next  observation  time  t^  ,  define 


T  s  *1+1  "  fci  ’ 


a  =  p  +  2td  +  t  P  • 
*XX  ^XU  LIU* 


(1) 

(2) 


and 


m  =  a  +  qx, 


x  =  xA  +  U^x  « 


e  =  zi+1  -  x. 


(3) 

(4) 

(5) 


and  denote  x(ti+1>  by  x  and  r^+1  by  r  for  brevity.  Then  it  follows  from 
standard  results  for  moments  of  Normal  random  variables  that,  given  z  ,...,z^, 


and 


E(e  )  sa+r+qx  , 
var(e2)  =  2  (m+r)2  +  3sx2, 

cov  (q,e  )  =  Sx. 


(6) 

(7) 

(8) 


By  assumption,  E(q/zQ, . . . .z^  =  q  (9) 

and 

var  (q/zQ,  • .  •  .Zj^)  =  s.  (10) 

If  at  this  point  the  (different)  approximation  is  made  of  treating  q  and 

2 

e  as  bivariate  Normal  random  variables  under  this  conditioning,  with  the  pre¬ 
ceding  mean  and  covariance  matrix  elements,  then  it  is  a  standard  result  that 

E(q/e2)  £  q  =  E(q)  +  COv(q^e--  [e2-E(e2)]  (11) 

var(e  ) 

and  2  2 

var(q/e2)  £  c  *  var(q)  -  -0V  1  ,  (12) 

var(e  ) 

where  the  conditioning  on  z  ,  ...,z  is  suppressed  in  the  notation.  Substi- 
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tuting  Eqs.  (6)  -  (10)  into  (11)  and  (12)  we  obtain 


q  =  q  + 


2(i2iL)  +  3s 


(t2  -(»«■>  _ .) 


and 


c  =  s  - 


2(— ) 


3s 


(13) 


(14) 


p 

Reducing  the  denominators  of  Eqs.  (13)  and  (14)  to  "2(m/f)  +  s"  would  convert 
these  equations  to  the  form  of  the  updating  step  in  the  standard  Kalman  filter 
for  constant  q  with 


and 


q  -  measurement  = 


2 

e 


-  (a+r) 

T 


2 

measurement  noise  variance  =  2(m/x) 


at  measurement  epoch  i+1.  This  change  is  justified  to  some  extent  by  the  fact 
that  m  really  depends  via  (2)  and  (3)  on  the  previous  estimates  of  q,  which 
are  used  to  compute  pxx>  Also,  the  resulting  Kalman  filter  estimate  of  q 
reduces  to  simply  the  average  of  the  quantities 

e  -  ( a-i-r ) 

T 

observed  at  each  updating  time  when  m,r  and  t  are  all  equal  at  each  step, 
which  is  essentially  the  former  estimation  procedure  for  q  described  in  Ref. 
[1],  and  is  consistent  with  the  justification  given  there  as  a  least-squares 
estimation  procedure.  For  unevenly  spaced  observations,  however,  this  new 
updating  procedure  has  the  desirable  property  of  estimating  q  as  a  weighted 
average  of  these  quantities,  the  weight  being  less  for  epochs  at  which  t  is 
3mall  and  the  variance  of  the  "q  -  measurement"  is  therefore  large.  Thus,  it 
gives  a  sort  of  weighted  least  -  squares  estimate  of  q  when  used  sequentially 
in  the  overall  ship  tracking  algorithm.  Rearranging  terms  gives  this  Kalman 
filter  update  step  as 
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A 


and 


A  -  St  2 

q  =  q  +  - ?  "  -  "y  [e  -  (ra+r)] 

2  m  +  st 


a  = 


2  m 


.  2  2 
2  m  +  St 


s. 


(15) 


(16) 


Since  it  doesn't  make  sense  to  estimate  maneuvering  intensity  with  fewer  than 
three  observations,  such  updating  is  started  on  the  third  observation  (for 
which  t  >o)  with  s  =  ®,  which  corresponds  to  using  q  =  o  and  replacing  the 
"filter  gain" 


_ ST _ 

2  2 
2«  +  ST* 

1  2 
by  —  in  Eq.  (15),  and  using  a  =  2  (ih/t)  for  Eq.  (16). 

A  standard  Kalman  filter  for  estimating  a  constant  q  would  use  the  value 
of  a  from  one  update  time  as  the  value  of  s  for  the  next  update.  Since  the 
approximate  value  for  q  can  really  change  over  extended  periods  of  time, 
however,  it  is  more  reasonable  to  postulate  a  "forgetting  time  constant"  T  and 
use  exponential  deweighting  between  update  times  by  computing  "s"  from  the 
preceding  'V  a3 


[^preceding  ’  2(tu,  -  tj  ]  *  *(tM  -  tj  ('7) 

A 

to  avoid  "locking  in"  to  a  value  of  q.  This  particular  deweighting  scheme 
corresponds  to  q  changing  between  the  current  and  preceding  observation  times 
according  to 


q  s  ~  §T  +  white  noise,  (18) 

where  the  white  noise  intensity  is  such  that  the  steady-state  variance  of  q  is 

2  r - - — r~,  the  value  it  would  have  if  initialized  at  the  current  update 

i+1  "  o 

time. 
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The  estimate  q  here  can  assume  negative  values,  whereas  q  must  by 
definition  be  nonnegative.  Herce  the  latest  value  of 

A 

max  (q,  o) 

is  always  used  as  the  maneuvering  intensity  parameter  in  the  other  parts  of 
the  overall  ship  tracking  algorithm. 

As  a  last  refinement,  if  the  quantity 

A 

a  =  max  {o,  max  (q,  o)  -  max  (q,  o) }  (19) 

is  positive,  it  is  assumed,  as  explained  in  Ref.  C31,  that  puu  (which  is 
computed  at  preceding  update  time)  should  have  been  larger  by  approximately 

•  ,  where  t  is  the  current  time.  This  would  increase  a  by 

-  o 

2  -2  2  2  _  _ 

—  (m  -  m  )  =  —  (m+m)(m-m) 

T  T 

if  it  were  being  initialized  at  this  point.  In  the  above  in  is  the  value  of  m 
which  would  result  from  (2)  and  (3)  if  Puu  were  increased  as  described,  i.e., 

2 

-  aT 

m  * m  +  t=t~  • 
o 

To  compensate  for  having  reduced  the  effects  of  the  correlation  between  q  and 
m  from  those  in  (13)  and  (14),  the  quantity  ot  is  also  added  to  m,  giving 

m  s  m  +  ot(1  +  )  , 

o 

which  would  increase  an  initialized  value  of  a  by 


So  =  2C2  m  +  at(1  +  ^-)]  ^(1+  —-)  . 

o  ”  o 

To  be  conservative  and  to  avoid  dividing  by  t,  we  therefore  add  the  quantity 

So  s  2[2  m  +  at]  7^-  (20) 

”  o 

to  the  result  of  Eq.  (16)  at  all  update  times  as  an  approximate  correction  for 
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having  ignored  in  the  development  of  (16)  some  of  the  interactions  between 
estimating  q  and  estimating  x  and  u  in  the  overall  ship  tracking  algorithm, 
namely  the  dependence  of  puu  on  q.  Of  course,  the  other  corrections  of  this 
sort  in  the  overall  algorithm,  which  are  for  parameters  of  the  conditional 
distribution  of  x  and  u  as  described  in  Ref.  C3l,  are  retained  as  before. 

4.  PROCEDURE  FOR  TWO-DIMENSIONAL  MOTION 

Only  planar  motion  in  rectangular  coordinates  is  discussed  here.  The 
algorithms  for  tracking  such  motion  can  easily  be  extended  to  tracking  on  a 
sphere,  as  described  in  Refs.  [1]  and  [2].  In  the  context  of  planar  motion, 


x  =  ship  position 

v  =  ship  velocity  (average),  and 

z =  observed  position  at  time  t^,  i  =  0,1,.... 


are  all  2-vectors  with  components  in  these  two  coordinates,  and  the  composite 
"motion  state"  vector  has  four  components.  The  conditional  distribution 

of  this  state  vector  is  approximated  by  the  tracking  algorithm  as  4-variate 
Normal,  whose  mean  and  covariance  matrix  are  denoted  here  in  bivariate 
partitions  as 


X(t) 

r 

a  nH 

V  ( t ) 

d  llU 

p:<t) ;  P3ct) 


for  generic  time  t. 


For  convenience,  we  also  denote  the  current  estimate  of  the  maneuvering 

intensity  parameter  as  Q  and  the  covariance  matrix  of  the  error  in  the  i-th 

observation  as  R^,  both  of  which  are  2x2  symmetric  positive  semi-definite 

matrices.  Also,  §.(i)  and  q  (i)  are  used  to  denote  interim  estimates  of  the 
d  c 

(scalar)  down-track  and  across-traek  maneuvering  intensities  created 
immediately  after  the  i-th  observation,  and  and  to  denote  corresponding 
variance  parameters,  in  accordance  with  the  notation  of  the  preceding  section. 


The  operation  of  the  overall  ship  tracking  algorithm,  with  the  improved 
adaptive  estimation  procedure  for  maneuvering,  can  be  summarized  as  follows 
for  the  case  in  which  the  reports  all  specify  a  position  which  is  localized  in 
both  dimensions: 


11 


Initialization 


Upon  the  receipt  of  the  initial  report  (zq,  Rq)  at  tine  tQ,  tracking  is 
started  with 


A 


x(t*)  a  Z 

0  o 

(2  -  vector) 

v*:>  *  »o 

(2x2  matrix) 

P2(tl}  *  0 

(2x2  matrix) 

P3(t? 


'  2  • 
x  ;  o 


x  a  user-specified  prior 
speed  variance 


q.(o)  =  q  (o)  a  0  (or  some  other  user-specified  positive 

value) 


From  time  t*  to  time  t.+,  ;  i  a  0,1,... 

_ i _ i+1 _ _ 

Tracking  proceeds  by  performing  the  summary  steps  of  Fig.  1  as  follows, 
where  "a"  denotes  a  replacement  operation  as  in  FORTRAN. 

Step  1  -  Propagate  state  distribution  to  t^  by  successively  computing: 


c  .  =  max  {q  .  (i) ,  0) 
a  a 


c  a  max  Iq  (i) ,  0} 
c  c 


0  a  tan 


-1/  V1  (v 

V2 


v.  and  V£  are  the  two  components 
of  v  (usually  local  east  and 
north  components  of  velocity 
estimate) 


q11 

s 

2 

c  .  cos  9  c 
d 

c 

2 

sin  9 

qi2 

s 

(0d  - 

c  )  sin© 
c 

cos© 

q22 

s 

2 

c  .  sin  e  c 
a 

c 

2 

COS  9 

'qn 

1  *12 

Q  = 

— 

4- - 

.q12 

1  q22_ 

(2x2  matrix) 
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'i+1 


X  =  X(t*)  +  V(t*)  T 

A  =  P^t*)  +  [P2(t^)  +  PgCt*)  ]t  +  P  (t*)T* 
M,  =  A  +  Qt 

M2  =  P2(ti)  +  P3(ti)r 


(2  -  vector) 


2x2  matrices 


«3  =  yy 


At  this  point,  the  conditional  distribution  of  the  state  at  time  t 
given  zq,...z ,  can  be  approximated  as  (4-variate)  Normal  with 


i+1  * 


and 


mean  = 


x 

v(t*) 


covariance  = 


and  ti+1  need  not  be  an  actual  observation  time  for  this  purpose,  but  may  be 
any  time  after  t^  for  which  the  projected  probability  distribution  is  desired. 


Step  2  -  Update  state  distribution  with  new  report: 

Let  z  and  R  denote  zJ  .  and  R.  .  henceforth.  Compute 

i+1  i+1  r 


x(t.+. ) 
i+I 


x  +  M1  (M1  +R)”1 (z  -  x) 


v(till)  s  v(ti)  +  M2  <M1  *  R)’1(z  - 


2-vectors 


P1(titl)  =  Mi  -  Mi(Mi  + 


P2(ti*1>  s  M2  "  M1(M1  +  R>”1m2  }  2x2  matrices 
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P3(ti*1)  =  M3  -  M2(M1  +  R)_1M2 


Step  3  -  Update  of  maneuvering  estimate. 


For  i=0, 1  only:  set  q.  (i+1)  =  qH  (i) 


and  q  (i+1)  =  q  (i). 
c  c 


Otherwise,  compute  in  turn: 


y  =  t  -  t 
i+1  o 


=  SlM.il  and 


L  cj 


=  aRn  ;  a  a 


f  cos  9  I  sin  9*1 
L-sin  9  <  cos  9j 


i 

cos  9  i  sin  9 
-sin  e  j  cos  9 


(z  -  x) 


a  .  (3)  =2  (b./x) 
d  d 


a  (3)  =  2  (b  /t) 
c  c 


q.  (3)  =  -  Ul  -  b,) 


2 


'd  d' 


%  (3)  ■  7  “o  -  b<=>  J 


for  i=2 


"\ 


s.  -  ,  ‘Ca.U)  -  2  (b./r)  1  »  2  (b  /t) 
dad  d 


T 

s  =  e  T[c  (i)  -  2  (b  /y)2]  +  2  (b  /y)2 
c  c  c  c 
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od(i+1 ) 


2  bd  »« 

2  bl  * 


a  (i+1) 
c 


2  *1  ». 

2  2 
2b  +  t  s„ 
c  c 


for  i  >_  3 
only, 

T  is  user- 
specified 
"forgetting  time 
constant". 


q.(i+1)  =  q.(i)  + 
a  a 


v 


5k2  +  .  2 
2bd  V 


-rd] 


q  (i+1)  =  q  (i) 
c  c 


V 


,h2  +  a  2 
2b  Sr 
c  c 


Cec  -  bc  -rc] 


«d  =  max  {0,  max  {qd(i+1 ) ,0}  -  max  (qd C i ) , 0 } } 

a  =  max  {0,  max  {q  (i+1),0)  -  max  {q  (i),0}} 
c  c  c 


o  (i+1)  a  a  (i+1)  +  2  —  [2  b  +  a  rl 
d  d  Y  d  d 

<*c 

a  (i+1)  -  a  (i+1)  +  2  —  [2  b  +  a  r3 
c  C  y  C  C 

Step  4  (performed  for  i>2  only)  -  Adjust  Kalman  filter  update  by  successively 
computing 


2  2 
d,.  =  a.  cos  9  +  a  sin  9 
lid  c 


d,_  =  (a.  -  a  )  sing  cos9 
12  d  c 


2  2 

d22  =  ad  Sin  9  +  ac  008  9 


d  '  d 
11  i  12 

*  —  *  +  “ 

L  d12  !  d22j 


(positive  semi-definite  2x2 
matrix) 
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(2x2  matrix) 


pjUj*,’  .  P3(ti:1) 


=  *  R(M1  *  8  ♦  ♦  R  .  DtT'Rt 

(2x2  matrix) 

v(ti^1)  =  v(ti*1)  +  i  D(M^  +  R  +  Dt) (z  -  x)  J 

I  2-vectors 

A  A  1  I 

x(t.*1)  =  x(t  +  R(M1  +Dt  +  R ) D ( M  j  +  R)  (z  -  x)t  ) 

The  adjustments  in  this  step  are  explained  in  Ref.  [3].  Corresponding 
adjustments  for  the  maneuvering  intensity  estimation  procedure  constitute  the 
last  two  computations  listed  here  under  Step  3. 

5 .  ALTERNATE  METHOD 

This  same  kind  of  improvement  can  also  be  applied  to  the  alternate  method 
of  Ref.  [3]  for  adaptively  estimating  maneuvering  intensity,  in  which  the 
estimated  maneuvering  is  contrained  to  be  statistically  isotropic,  but  which 
allows  the  inclusion  of  bearing  -  only  input  reports  in  the  resulting 
tracking  algorithm.  In  this  case,  only  a  single  (non-negative)  scalar 
maneuvering  intensity  parameter  q  is  estimated,  and  the  overall  tracking 
algorithm  can  be  summarized  as  follows  for  position  reports: 

Initialization  -  as  in  the  algorithm  of  the  preceding  section,  except  that 
only  a  single  value  q(o)  is  specified,  instead  of  qd(o)  and  qc(o). 


From  t*  to  t^  ;  i  =  0,  1,  ...: 


Step  1  -  calculate  the  2x2  matrix  Q  by  the  steps 


a 

q  =  max  {q(i),  0} 
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Then  continue  as  in  the  preceding  section. 


Step  2  -  as  in  the  preceding  section. 


Step  3  -  For  i=0,1,  just  set  q(i+1)  =  q(i). 
Otherwise,  compute  in  turn: 


sn :  gi2 

_g12  !  g22  J 


=  A  +  R 


(2x2  matrix) 


g11  +  g22 


(g11  “  g22)  *  4  8 
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a2  =  1  (a2  *  b2) 


.2  2  K2 

b  =  a  -  b 


cos  8  = 


sin  s  = 


M2 


/(a*  -  gn>*  ♦  g,2 


Note:  a  and  b  are  the  semi-major  and  semi-minor  axes  of  the 
the  bivariate  Normal  distribution  with  covariance  matrix  A  + 
between  the  semi-major  axis  and  the  1-coordinate  axis. 


ea 

'  coss  |  sins  " 

?v 

-sins'  r  cos~s” 

'  J 

(z-x) 


(2-vector) 


sigma-ellipse  for 
R;  8  is  the  angle 


(21) 
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2  2 
h  =  g^  sin  s  -  2g10  sin  8  cos  6  +  g,5  cos  8  +  qr 


M2 


s22 


n  =  b  +  qx 


k  =  a  +  qx 


a(3)  =  2(-)2 

T 


q(3)  =  -(.cl  -  n) 
x  b 


for  i=2  only 


(22) 


s  =  e  T  ta(i)  -  2(-)2]  +  2(-)2 
y  y 


c(i+1)  = 


2h2s 


2h  +  st 


S  for  i>3  only 


q(i+1)  =  q(i)  +  — 2^ - 2  ^eb  “ 

2h  +  sx  J 


o(i+1)  = 


2  k 


2k  +  sx 


r(i+1) 


(23) 


q( i+1 )  =  q(i+1 )  +  — 12 - -  k] 

2k  +  sx 


(24) 


a  s  max  {0,  max  {q(i+1),0}  -  max  {q(i),0}} 


a(i+1 )  2  a(i+1 )  +  2  —  [2h  ♦  axl 


Step  4  -  calculate  the  2x2  matrix  D  as 
then  continue  as  in  the  preceding  section. 


•  ■  m  ■ 


The  advantage  of  this  alternate  method  is  that  it  can  make  effective  use 
of  bearing-only  reports  of  ship  location  (one  of  which  specifies  a  time,  an 
observed  bearing  from  a  specified  point  at  that  time,  and  a  bearing  variance). 
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I 


k 


or  what  is  almost  the  same  thing,  position  reports  with  long,  narrow  error 

ellipses.  One  simple  way  of  doing  this  with  such  a  bearing  report  is  to 

replace  it  with  an  approximately  equivalent  position  report,  which  would  have 

a  large  variance  component  along  the  line-of-bearing,  and  use  the  procedure  as 

just  described.  The  best  way  of  constructing  an  equivalent  position  report 

depends  on  the  bearing  sensor  range,  how  much  the  projected  ship  position 

distribution  at  the  time  of  the  report  "overlaps"  the  bearing  sensor  location, 

and  the  distance  of  the  estimated  ship  position  at  that  time  from  the  3ensor 

location.  The  main  idea  is  to  have  an  appropriate  containment  ellipse  for 

this  (bivariate  Normally  distributed)  position  report,  say  the  two-sigma  (86%) 

ellipse,  match  the  corresponding  wedge-shaped  confidence  region  of  the  bearing 

report  as  well  as  possible  within  the  range  of  likely  projected  ship 

locations.  As  a  numerically  more  efficient  alternative  when  the  predicted 

ship  range  from  the  bearing  sensor  is  large  compared  to  the  uncertainty  in 

this  prediction,  one  could  also  use  a  bearing-only  report  as  a  scalar  measure- 

2  2 

ment  of  the  cross-bearing  component  of  ship  position  with  variance  r  5  ,  where 

2 

r  is  the  projected  ship  range  (from  the  sensor)  and  6  is  the  bearing  variance 

in  radians.  This  would  require  the  replacement  of  Step  2  (in  this  section)  by 

the  appropriate  Kalman  filter  updating  procedure  for  such  a  single-component 

2  2  2 

measurement.  It  is  also  necessary  to  alter  Step  3  by  computing  b  as  r  6 
plus  the  cross-bearing  component  of  computing  as  the  cross-bearing 
component  of  the  observed  minus  predicted  ship  position  (as  evaluated  from  the 
observed  bearing  at  the  predicted  ship  range)  in  lieu  of  Eq.  (21),  and  skip¬ 
ping  the  computations  of  Eqs.  ( 22) —(24 )  and  those  prior  to  Eq.  (21). 

For  position  reports  which  are  well-localized  in  both  dimensions,  the  use 
of  this  alternate  tracking  algorithm  would  probably  not  cause  much  loss  of 
precision  over  the  algorithm  of  the  preceding  section.  It  was  originally 
thought  that  allowing  the  estimates  of  the  down-track  and  cross-track 
components  of  the  maneuvering  intensiity  to  be  different,  as  they  may  be  in 
the  algorithm  of  the  preceding  section,  would  often  lead  to  more  precise 
estimates  of  ship  motion  when  a  ship  zigzags  about  an  average  course,  or 
changes  speed  a  lot  but  stay 3  on  the  same  heading.  In  practice,  however,  the 
estimated  maneuvering  is  usually  fairly  isotropic  anyway. 


I 
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6.  DISCUSSION 


To  summarize  the  basic  usage  of  the  two  ship  tracking  algorithm  variants 

presented  here  for  the  case  of  position  reports,  the  user  specifies  the  two 
(scalar)  parameters: 

X  =  Prior  estimate  of  ship  speed  (i.e.,  average  speed  for  anticipated 
target  population).  Making  x  too  large  doesn't  matter  much  except 
when  this  algorithm  is  used  as  part  of  a  multitarget  tracker,  but 
making  x  too  small  does. 

T  =  Average  length  of  time  over  which  a  ship's  maneuvering  behavior  is 

judged  to  remain  statistically  the  same.  This  value  can  often  be  made 
infinite  without  serious  consequences;  another  reasonable  possibility 
is  to  set  T  =  /A/x,  where  A  is  the  area  of  the  surveillance  region  in 
which  ships  are  being  tracked. 

Also,  it  is  sometimes  slightly  advantageous  to  specify  .iome  nonzero 

A  a  A 

initial  value  for  qd(0)  and  qc(0),  or  q(0)  in  the  case  of  the  variant  of 
Section  5.  For  either  variant,  the  algorithm  then  operates  on  input 
consisting  of  a  time-ordered  sequence  of  ship  position  reports,  each  of  which 
contains  the  following  six  data: 

1.  time  of  report 

2-3.  two  coordinates  of  observed  ship  position  at  that  time 

4-6.  the  three  independent  components  of  the  (symmetric  2x2)  covariance 
matrix  for  the  error  in  this  observed  position. 

As  output,  it  provides  (a)  estimates  of  ship  position  and  average  velocity,  at 
present  or  future  times,  (b)  error  covariance  matrices  (4x4)  that  correspond 
to  these  estimates,  and  (c)  estimates  and  corresponding  error  covariances  of 
maneuverability  parameters  (two  in  the  variant  of  Section  4,  one  in  the 
variant  of  Section  5). 

Normally,  the  only  outputs  of  interest  are  the  position  components  of  the 
estimated  position-velocity  state  vector,  and  the  corresponding  2x2  error 
covariance  matrix.  It  is  also  usual  to  regard  these  errors  and  the  position 
observation  errors  in  the  input  reports  as  bivariate  Normally  distributed,  and 
to  use  equivalent  parameters  specifying  the  two-sigma  ellipses  (86%  contain¬ 
ment  ellipses)  of  these  distributions  in  place  of  their  covariance  matrices. 
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The  ellipse  parameters  normally  used  for  this  purpose  are  the  semi-major  and 
semi-minor  axis  lengths  and  the  angle  of  the  major  axis  measured  clockwise 
from  local  north.  As  examples  of  what  this  all  amounts  to  in  practice, 
experimental  FORTRAN  implementations  of  both  variants  of  this  overall  ship 
tracking  algorithm  are  listed  in  the  appendix  for  the  usual  case  in  which  the 
input  and  output  error  distributions  are  specified  in  terms  of  containment 
ellipse  parameters  instead  of  covariance  matrix  components.  The  input  and 
output  are  from  a  terminal  in  these  implementations. 

Figures  2  and  3  show  an  example  of  the  comparative  performances  of  this 
tracking  algorithm  (the  version  of  Section  4  here)  and  the  corresponding  one 
of  Ref.  [3],  which  uses  the  former  method  of  adaptive  maneuvering  estimation. 
The  parameter  values  q  .(0)  =  q  (0)  =  0  and  T  s  »  were  used  in  each  case.  For 
clarity  in  these  figures,  the  two-sigma  (.86%  containment)  ellipses  for  the 
output  location  estimates  are  displayed  only  for  selected  times.  The  new 
algorithm  is  not  shown  to  its  best  advantage  in  this  example  because  the 
position  reports  are  all  evenly  spaced  in  time.  Even  so,  however,  the 
accuracy  of  the  estimated  ship  positions  is  as  good  as  that  of  the  former 
algorithm,  and  the  containment  ellipses  generated  by  the  new  algorithm  are 
generally  somewhat  larger,  which  is  an  indication  of  better  statistical 
consistency  in  light  of  the  computational  experience  with  actual  tracking  data 
mentioned  earlier. 

The  variant  of  Section  5  here  has  also  been  embedded  in  a  multi-target 
tracker,  in  the  manner  described  in  Ref.  [4],  as  a  replacement  for  the 
corresponding  single-target  tracking  method  of  Ref.  [33,  which  estimates  only 
a  scalar  intensity  for  statistically  isotropic  maneuvering.  The  input  data  in 
this  application  were  mostly  bearing-only  reports  which  occurred  in  sporadic 
bursts  for  each  ship,  and  thus  were  very  unevenly  spaced  in  time.  In  this 
case  there  was  a  substantial  improvement  in  performance  with  the  tracker 
incorporating  the  new  adaptive  maneuvering  estimation  procedure  described  in 
this  report. 


Fig.  2  —  Tracking  performance  with  former  procedure 


Fig.  3  —  Tracking  performance  with  procedure  of  Section  4 
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APPENDIX 


PROGRAM  LISTINGS 


C  KALMAN  FILTER  with  a  J  APT  IV£  DRIVING  NOISE 

C 

C  FOR  TRACKING  IN  RECTANGULAR  I ,  Y-COORO INATES  wITM  POSITION 
C  REPORTS  -M ICW  are  MEaninGFuLl"  LOCALIZED  IN  TWO  OINEnSIOnS, 
C 

c 

C  INITIALIZATION 

c 

c 

c  enter  user-specifieo  parameter  values  fron  terminal 
c 

ACCEPT"PRIOR  expected  ship  speeo  »  ",a 

c 

C  COMPUTE  INITIAL  VALUES 

c 

Ni«a. 

WC»«. 

ci»b, 

CC»R . 

u«e. 

VaR, 

PXU*0. 

P XVaR. 

PVU»R. 

PVVafl, 

Puu«,a«ii*a 

PUV"<J, 

PVVaPUU 

9XXa0, 

OXYap, 

OYY»0, 

C 

C 

C  9EGIN  TRACKING 

C 

C 

c  enter  ano  process  initial  report 
c 

C  IR  »  INDICATOR!  ZERO  IF  "REPORT"  JUST  SPECIFIES  A  TIME  AT 

C  -M ICM  OEAO'RECKONED  output  IS  OESIRED. 

C  POSITIVE  if  IT  is  AN  actual  report  of 

C  OBSERVED  POSITION,  NEGATIVE  TO 

C  TERMINATE  TRACKING, 

c 

C  FOR  POSITIVE  IR I 

C  T  a  TIME, 

C  ZX.ZV  a  uaSERVEO  POSITION  IN  X , Y-COORO IN AYES , 

C  SNA  a  SENiNAJOR  axis  OF  B8X  CONTAINMENT  ELLIPSE  FOR 
C  OBSERVATION, 

C  SMI  a  SEMIMINOR  axis  OF  CONTAINMENT  ELLIPSE, 

C  TMT  a  ORIENTATION  OF  MAJOR  AXIS  (DEGREES  CLOCKalSE  FROM 
C  Y-AXIS). 

c 

ACCEPTOR,  T, ZX.ZV, SNA, SMI, TNT  a  ",  IR,  T,  ZX,  ZY,  SMA,  SM I ,  TMT 
IFCIR.LE.R)  GU  TO  99 
I  *R 

TmTaTmT/SF, J 
xaZ( 

Y  aZY 

STaSlN(TMT) 

CTaCOS(TMT) 

SM  A  aS*"  A  *SM  A 
SMlaJMtaSMl 
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PXXa,25a(SPAaST«ST*S0I*CT.CT) 

PXYa,2SaSTaCTa(SMA-SHIl 

PYYa.25a{SRAaCT.CT«SHI»$TaST) 

TNaT 

tl«t 

TLU«T 
GU  TO  40 
C 
C 

C  CONTINUE  TRACKING 

C 

c 

c  enter  next  report 
c 

90  ACCEPTOR, T,ZX,ZY,SHA,SMI,ThT  •  " , IR , T , ZX , Z Y , SHA , SHI , TMT 
IP  tIR.LT.RJ  GO  TO  99 
IPCIR.EQ.0J  GO  TO  10 
ThT»ThT/S7,J 
STaSInCTNTJ 
CT«COSCTHT) 

$ha«Sha*sha 

ShI»SnI«Sh1 

RxX»,25»(SHA»ST.ST*SMI.CT«CTJ 

RXY».25«ST«CWSRA-SHI) 

RTY».25*tSHA*CT*CT*SHI*ST*STJ 

c 

C  PROJECT  STATE  OENSITY  TO  CURRENT  TIRE 
C 

10  TAUaT-TL 

tl*t 

Xo*X*U«TAU 

yb«y*v«tau 

5XXiPXX*(2,«PXU*0XX)»TAu»PUU«TAU»TAU 
GXY.PXY*(PXV»PYU»QXY) aTAU*PUVaTAU«TAU 
GYY«PYY*(2,»PYV*QYY)«TAU*PVV«TAU*TAU 
gxu>pxu*puu«Tau 

GXVaPXVaPuVaTAU 

gyu«pyu*puv«tau 

GYV«PYV*PyV«TAU 
IPCIR.NE.0J  GO  TO  S0 
X»X0 
Y»Y0 
PXXaGXX 
PxYaGXY 
PYY»GYY 
PXUaGXU 
PXVaGXV 
PYUaGYU 
PYYaGYV 
GO  TO  40 
C 

C  UPOATE  STATE  OEnSITY  »ITH  N£R  REPORT 
C 

60  TAUaT-TLU 
TLUaT 
ExaZX-Xo 
EYaZY-YB 
Gt "GXXaRXX 
G2*GXY*RXY 
G3»GYY*ryy 
OETaGl»G3-G2«G2 
OETIa«, 

IPCOET.GT.0.J  OETIal./OET 
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«XX<G3<0fcTI 

HXT<<U2<0ETI 

HVVlGl<OETI 

X<XB<(GXX<HXX<GXY<MXY3*EX<(GXX<MXY<GXY*hYYJ*EY 
Viva<CGXY<HXX<GYr <HXV3<EX<CGXY <H XT  <GYY<HVTJ<EY 
U <U<CUXU<HXX <GYU<HXYJ<EX<CGXU<HXY<GYU<HYYJ<EY 
V<V<CUXV<HXX<GYV<HXY3<EX< CGXV<HXY<GYV<HYYJ<EY 
PxX<GXX-GXX<GXX<HXX-2, <GXX<GXY<HXY<GXY<UXY<HYY 
PXY<GXY»GXX<GXY<HXX<(GXX<GYY<GXY*GXY3*HXY<GXY<GYY<NYY 
PYY<GYY-GxY<GXV<HXX"2.<GXY*GYY<HXY<GYY*UYY<hYY 

pxu<gxu-gxx<gxu<mxx<(gxx<gyu*gxy<gxuJ <hxy-gxy<gyu<myy 

PXV<GXV<G1X<GXV<HXX<CGXX<GYV<GVY<GXV3<MXY <GXY<GYV<HYY 
PYU«GYU-G XY<GXU<HXX<IGX Y<GYU<6YY<GXU3  <HXY<GYY<GYU<HYY 
PYV«GYV-GXY<GXV<HXX<(GXY <GYY<GYY<GXV3<HXY<GYY<GYV<HYY 
PUU<PUU<GXU<GXU<H XX <2, <GXU<GYU<HX Y<GYU<GYU<HYY 
PUV <PUV<GXU<GXV<HXX<CGX U<GYV<GYU<GXVJ  <HXY<GYU<GYV<MYY 
PW  <PVV<G  XV  <GXV<HXX<2,  <GXV<GYV<HXY  <GYV<GYV<HYY 
C 

C  UP04TE  USIYJHG  NOISE  ESTlHXTE  4N0  40JUST  ESTIH4T0N  P4H4HETEHS 
C 

IP (TXU.UE.P.3  GO  TO  42 
Tt«T<lN 

IP  CTE.lt. a. 3  SO  Tu  43 
I<I*1 

IFCI.tQ.13  GO  TO  43 

U2«u<g 

V2 <Y*V 

3U*u2<V2 

Q1NV<U. 

IF (SQ.GT .3.3  UINV<1,/SQ 

U2<U2<<JINV 

V2<V2<<3INV 

UV<U<V<UINV 

SIT<QINV<CU<EX<V<EY3<<2 
SCT<U1NV<IU<EY<V<EXJ<<2 
GIT<U2<GXX <2. <UV<GXY<V2<6YY 
GCT<V2<GXX -2, <UV<GXY<U2<GYY 
0IT<U2<Gl<2. <UV<G2<V2<63 
8CT<V2<Gl<2. <UV<G2<U2<G3 
IFCI.UT.aj  GO  TO  33 
GNI<1,/T<U 
GNCsGNl 

VHI<2, <GIT<GIT<GNI<GNI 
VHC <2, <GCT<GCT<GNC<GNC 
GO  TO  33 

33  0EN<2, <GIT<GIT<VRI<T4U<T4U 
OtNI<*. 

IF COEn.GT.3.3  OENl«l,/OEN 

GNI<VHI<TXU<OENI 

VKI<VMJ<(1, <VWI<T4U<T4U<0EN13 

06N»2, <GCT<GCT<VSC<T4U<T  XU 

OENI<u, 

I F  CUE N  , GT . 2 , J  D£NI<1./0EN 
GNC<VHC<TXU<OENI 
VHC<VHC<C1, <VHC<TXU<T4U<0ENIJ 
33  C1<CI<GNI<CSIT<8IT3 
CC<CC<GNC<CSCT<8CT3 

0  I  <  H  I 

OC«HC 

Nl<8. 

HC<3. 

IFCCI.GT.8.3  h I <C 1 
I< CCC.GT.3, J  HC»CC 
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cc»«c-oc 

QXX>U2*hI*V2*mC 

qxy*uv* 

0YY*V2.hI«U2«*C 

PI««. 

PC»4. 

IFCOJ.GT..).)  P!«OI/TE 
JFCDC.GT.U.)  PC*OC/TE 
VKI«yi<:*2..p:«  (PI.Tfe«T*U*2.»GITJ 
VhC»V>«C*2.*PC«  tPC*Tt*TAU*2.*GCT) 

C 

C  ADJUST  PANAMETEKS  OF  state  DENSITY 

c 

EZ*"Xx*tX»HXY»£Y 

t«*Mxr»tx*^TT«eT 

1XX»U2«PI»V2»AC 

HXY*UV«(PI-PC) 

HYYaV2*PI*U2*PC 
PUU«PUU**XX 
PuV*PUV*HY 
PWaPyVANYY 
A  X  X  AH  1 X  *  T  AU 

axy*«iy«tau 

ayy»«yy«tau 

GXX«A1X»TE 

Gxy«aiy«te 

GYYaAYY#T£ 

Gl*ul*Gxx 

G2*G2*GXY 

G3«G3*GYY 

DtT«Gl»G3-G2*G2 

0ETI*4. 

IPtDET.GT.a.)  OETI«l,/DET 
hxx*G3*0ETI 
MXY«-G2«0ETI 
hyY«G1*OETI 

U»U» fAXx»«XX*AXY«HXY) •EX*IAXX»MXY»AXY«NYY J«EY 

V«V»(AXY«nXX*AYY*HXY)*EX»CAXY*HXV»ATT*MYY) *EY 

Gl*HXX«*XX*rtlY«MXY 

G2*PXX»mXY*SXY«MYY 

G3*i»XY»«XI*4YY*HXY 

G4*NXY»NIY*HYY«HYY 

PX<>PXX*GXX»G1*G1*2.*GXY*G1*G2*GYY*G2*G2 
PXYaPXY*GXX*Gl*G3*GlY*CGi*G4»G2*G3)*GYY*G2*C  1 
PYY>PYY*GXX*G3*G3*2a*GXY*G3«G4»GYY*G4*G4 
X*X*GXX«G1*EZ*GXY*(G1«E**G2*EZ)*GYY»G2*£m 
Y«Y*Gxx*G3*£Z*GXY*(G3*Ei>*G4*EZ}*GYY*G4*E4 


OUTPUT 


4V  C 1 aPXXaPYY 

C2*SGNT{CPXX«PYY)*«2*4,*PXY«PXY} 

Cl*. 3* (C1*C23 

C2*C1«C2 

S.«A*2,»SOHTtClJ 

SM»2.»30rfTCC2) 

IF(PXT,NE,t».)  GO  TO  7Z 
THT*fl, 
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IPtPXX.LT.PYY)  TMT«90. 

50  TO  Sit 

79  TMT«S7,3»ATAN((PXX«Cn/PXY)*9t!, 

OISPLAY  OUTPUT  ON  TERMINAL 

IR  i  AdOVE-OESCRIBED  IN0ICAT0R  POP  "REPORT*  JUST  PROCESSEO, 
T  *  CURRENT  TIRE  (»  TIME  OF  This  REPORTS, 

X,Y  ■  CURRENT  ESTIMATED  POSITION  IN  X  ,  Y-CQORUI NA TES , 

5" A  ■  SERI  MAJOR  AXIS  OP  B6X  CONTAINMENT  ELLIPSE  F0« 

CURRENT  POSITION  ESTIMATE, 

SRI  •  SfcMlMINOR  AXIS  OP  CONTAINMENT  ELLIPSE. 

TMT  ■  ORIENTATION  OP  MAJOR  AXIS  (DECREES  CL0C*«I5E 
PROM  Y.AXIS, 

40  TYPE"  " 

TYPE"IR  •  ",I R 
TYPE'T  <  »,T 
TYPt"X,Y  •  * , X , Y 
TYP£"SMA,SMI  »  ",  SNA, SMI 
TYPt’TMT  ■  ",TMT 
TYPE*  * 

C 

C 

C  PROCESS  NEXT  REPORT  (OR  OEAO-RECaON  TIMES,  IP  ANY 

C 

c 

SO  TO  94 
C 
C 

C  TERMINATION 

C 

C 

99  CONTINUE 
ENO 
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MI.HAN  FILTER  nITM  AUAPTIVt  DRIVING  NOISE 

FOR  TRACKING  In  RECTANGULAR  X,  Y-COORDInaTES  »ITh  POSITION 
REPORTS  <mICn  ARE  HEANINGFULLT  LOCALIZED  In  ONLV  ONE 
OIMEnSION,  Such  AS  APPROXIMATIONS  OF  BEAR ING-ONL  V  REPORTS, 


initialization 


enter  user-specified  parameter  values  from  terminal 

ACCtPT»PRIOR  EXPECTED  SmIP  SPEED  •  ".B 

c 

c  COMPUTE  INITIAL  values 
C 

CI*v . 

U»P  , 

V»t6  . 

PxU»0. 

PXVatt, 

PTU«0, 

prv»n, 

puu«.a«o»B 

PuVaM, 

PVVaPLU 
QaH  , 

c 
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c 

C  6EGIN  TKAOING 

C 

C 

C  EnTeR  and  PROCESS  IMT1AL  REPORT 
C 

C  IR  ■  INDICATOR!  ZERO  IF  "REPORT"  JUST  SPECIFIES  A  TIME  AT 

C  kHICh  OEAO-MECKOnED  OUTPUT  IS  OESTREO. 

C  POSITIVE  IF  IT  IS  AN  ACTUAL  REPORT  OP 

C  OBSERVED  POSITION,  NEGATIVE  TO 

C  TERMINATE  TRACKING, 

C 

C  PQR  POSITIVE  IRl 

C  T  •  TIME. 

C  ZX.ZY  *  OBSERVED  POSITION  IN  X , Y-COORO INATtS , 

C  SNA  ■  SEmImaJOr  AXIS  Of  8SX  CONTAINMENT  ELLIPSE  FOR 
C  OBSERVATION, 

C  SMI  *  SEmImInOr  AXIS  OP  CONTAINMENT  ELLIPSE. 

C  TNT  •  URIEnTaTIOn  OF  MAJOR  Alls  [DEGREES  CLOCKrISE  FROM 
C  T-AXIS). 

C 

ACCEPTOR, T.Zx.ZY, SNA, SMt.THT  *  *  ,  IR  ,  T  ,  ZX  ,  ZT  ,  SNA ,  SMI ,  TMT 
IFC1H.LE.U]  GO  TO  99 
I  PM 

ThT  pThT/ 57 , 3 

x»Z» 

Y  •  Z  Y 

ST>SlN(TMT) 

CTpCOSCTmTI 
Snap Sma* Sma 
SMl«SMl«SMl 

PXXp.25* (SMA»ST»ST*SMI«CT«CT) 

PXYp.2S«ST«CTp<SMa-SMI3 

PYYp,25«CSMa»CT«CT*SMI*STp$TJ 

TnpT 

Tl«T 

TlUpT 

GO  10  4k> 

C 

C 

C  CONTINUE  TRACKING 

C 

C 

C  ENTtR  NEXT  report 

c 

9U  ACCEPTOR, T.Zx.ZY, 5mA, SnI.THT  p  ",IH,T, ZX.ZT, SNA, SMJ.TMT 
IPCIR.lt. BJ  GO  TO  99 
IF (IR.Eu.u)  GO  TO  IV 
TmTpTmT/37,4 
STpSInCTMTI 
CTpCOSCTmTJ 
SmapSmapSma 
5m  I  »SM  I  »SmI 

RXX»,25*CSMA"ST*STpSMI*CT*CT) 

RXYp,25*ST"CT»(SMA-SMI} 

RY*p.<!5"CSMa»CT.CT*SMI.ST«ST) 

C 

C  PROJECT  STATE  DENSITY  TO  CURRENT  TIME 

c 

\e  taupt-tl 

TLPT 

XbPl*U»TAU 

YBpY»v«TAu 


uu  U  OOU 


GA*FXX*2. •PXU»TAu*PUU*TaU*TAU 
Gb«PYY*2,«PYV*TAU*PVV*TAU*TA0 
GXX  *GA*G*  r  iu 

GxY*PXY*(PXv*PYU)*TAU*PuV*TAU*TAU 

4YY*Go*U»  FAu 

GXU*PXU*PUO»TAli 

6xv*pxv*Puv*Tau 

GYU*PYU*PUV*TAU 

GYV«PYV*PW.Tau 

IP  GO  TO  60 

XsXo 

Y»YB 

PXXsGXX 

PXY»GXY 

P»Y»GYY 

PXUsGXU 

PxVbGXV 

PYU*GYO 

P»Y»GYV 

GO  TO  AO 

UPDATE  STATE  OEnSITY  »1Th  NE*  REPORT 

60  TAO»T«TLU 
TLLiaT 
Ex«lx-xe 
£v*2*-Yo 
G1*GXX«hXx 
G2«GXY*«XY 
GJ«uYY*RYY 
0ET«G1»G3-G2*G2 
D£TI«n, 

IFCUET.GT.0,3  0£TI«1  ,/OET 
MXX*G3'»uETI 
HAY«-G2*0ETI 
hYY»G1*uETI 

X»X6«(GxX.HXX*GXY*HXY)*EX*CGXX*MXY*GXY*MVY)*EY 
Y*Yo*lGXV *HXX*GYY.HXY) •£X*CGXY*hXY*GYY*HYY)*EY 
U«U*(GXu*"*X*GYu«hXYJ •EX*(GXU*MXY*GYU*MYY)*EV 
v«v»(U*v»MXX*GYV»NXY)»£X*(GXV*MXY*GYV»MYYJ  *£  Y 
PXX*GXX-GXX*GXX.HXX-2.«GXX»GXY*MXY-GXY*GXY*MYY 
PXY»GxY-GXX«GXY.MXX-(GXX*GYY*GXY*GXY> *mxy-Gxy*gyy*myy 
PYV»GYV-GXY*GXV«MXX-a,.GXY*GYY*MXY-GYY*GYY*MYY 
PXU«GXU-GXX«GXU*HXX- tGXX«GYU*GXY*GXU) *mx y-gx y *g yu«m y y 
PXVBGXV -G XX. GXV*MXX-(GXX*GYV*GXY*GXV] *MXY-GXY*GYV«MYY 
PYU*GYU-GXY*Gau*MXX-(GXY»GVU*GYY*GXU)*HXY-GYY*GYU*MYY 
PYV*GYV-GXY*GXV*HXX-CGXY*GYV*GYY*GXVJ •MXY*GYY*GYV*MYY 
PUU*PUU-GXU*GXU*M XX -J ,»GXU*GYU»HXY-GYU»GYU«HYY 
puv*puv-Gxu»Gxv«mxx-(gxu*gvv*g.yu*Gxv)*mxy-gvu*GYV*myy 
PVV*PVV«GXV*GXV»MXX-2,*GXV*GYV*MXY-GYV*GYV*HYY 

UPDATE  DRIVING  NOISE  ESTIMATE  AND  ADJUST  ESTIMATOR  PARAMETERS 

IF CTau.uE.B.J  GO  TO  40 
TE»Y-TN 

IFITE.Lfc.0.)  GO  TO  40 
1*1*1 

XFCI.EQ.l)  GO  TO  40 
01*G 

G1*GA*RXX 
G5*GB*RyY 
C 1  *G 1 *G3 

C2*SG*TI  tGl-G31**2*4..G2*G21 
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non 


C1«.5«(C1*C2J 

C2«C1-C< 

Cb*Cl-Gl 

0tN*Cb*C8»G2*t»2 

ObNlit.. 

IK(0En,uT.V.)  0£NI«I./0£n 
IF (C2.GT,.90*Ct)  C8«l. 

IF0.2.GT..9**C1J  o£NI»l. 
EA»0Ef<I*{Cr«£»»G2»EYJ*«'<i 
£b»U£Nl*(0»£Y-G2»EX  J**2 

Hi(G2*G2«GXX»2,«G2*CB«G*Y*Cb»CB»GYY)»UENI 

hh*C2*0«  T  40 

IF  (I. GY. 2)  GO  TO  3P 

GNl«l ,FTAU 
Vnl«2,»h«"«GNI«GNl 
SO  TO  Sf 

3^  0c*»2.»n*"«'Yf<I«T»U«TAU 

IF (U£n ,GT .9.5  OEhlal ./OtN 

GNl*v«I«TAO»OtMI 

VxI»V«I«U.-*p<1«Tau»TAU«DENI) 

St*  CI*CI*GNl«(tS-H«3 
Q«F. 

IFICI.GT.il.)  <J»CI 
ha»C1*0»TAU 

0£n«2,»ha»pa»^#I»TAU»TAU 

OcNIan*. 

IF  tuEi»,GT.^.  J  0£KI»1  ,/DEN 
GkI*vHI«Tau*OENI 
VRI»VmI«(1.*>VWI»TAU»TAU«0ENI) 
CI«CI*G«l»(fcA-HAI 
0*0. 

IFtCI.GT.M.)  o«ci 

DI»Q-OI 

Pl«ie. 

IFlOI.GT.b.)  PI«OI/TE 
VHI»v«I*2,«Pl» (PI«TE*TAu*2,«*0 

AOJUST  PARAMETERS  OF  STATt  DENSITY 

fcZ«N«*aE**FiAY»EY 

E*«ri*Y»tX*NYY«EY 

Pull»Puu*PI 

PW»P*V»PI 

a»PI«Tau 

G*A»Tt 

Gl*bt*G 

G3*G3«G 

0£T»G1*G3-G2*G2 

DETI«U. 

IFC0ET.GT.9.J  OET I  *  1 .  /OET 

HXXsG3*0£Tl 

hxy«.G2*0ET1 

HYYaGl *OETI 

U*U*A*(nXX«£«»HXY*EY] 

V«V»A«(nXYa£X*HYY»EY) 

Gl«RXX»»1XX»NXY«HXY 
Gi»KXX«NXY*»XY«MYY 
G3»RXY»m»i»«YY«MXY 
G««*>XY«MXY*hYY«HYY 
Pll>PXX*G* ( G 1 *G 1 *G2  »G2) 
plY«PXY*G« (GI *G3»G2»G«) 

PYY»P»Y*G« (G3*G3*G4»G«) 
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X.»*G.(G1«£Z»&2«E<0 

yay.GaCGS«E2*G4«Ea) 

C 

C 

C  Output 

C 

C 

40  Cl»PXX*PYY 

C2*SQRT(  (PXX-PVY) ••2*4, »PXT«PXT) 

C  l  ■  ,S* (C1*C2) 

C2-CI-C2 
SP4«2,«aOPTlCn 
S* I  »2 .  «sQk T  (C2) 

IT  (PXY  nE.0.)  GO  TO  7(9 
ThTaS), 

IKPY4.GT.PYY)  THTa90, 

Gb  TO  a* 

70  T«T«S7. J«ATAN((PXX-Cn/PXT)*90, 

C 

C  DISPLAY  OUTPUT  QN  TERMINAL 
C 

C  I«  «  AaOvE-PtSCRISED  InOICATOR  FOR  "REPORT"  JUST  PROCESSED. 
C  T  ■  CURRtNT  TIME  (>  TIME  OF  THIS  REPORT) , 

C  X i Y  •  CURRENT  t5T I  MATEO  POSITION  IN  X , Y-COURO INATES . 

C  SNA  «  SEMnAJOR  AXIS  OF  06X  CONTAINMENT  ELLIPSE  FOR 
C  CU«»tNT  POSITION  ESTIMATE. 

C  SMI  a  SEmImIrOR  AXIS  OF  CONTAINMENT  ELLIPSE, 

C  TMT  a  uR IEnT AT ICn  OF  MAJOR  AXIS  (DEGREES  CLQCKaISE 
C  FROM  » •  A,x I S , 

C 

00  TYPE"  • 

TYPt"I«  •  ",IR 
TYPt"T  a  ",T 
TfPfX.Y  a  ",x.Y 
TYPt’SMA.SMI  A  ".SMA.SMI 
TYPE’TMT  a  »#TMT 
TYPE"  " 

C 

C 

C  PROCtSS  NEXT  HEPOHT  (OR  OEAO-»ECxON  TIME),  IF  ANY 

C 

C 

Gb  TO  90 
C 

c 

C  TERMINATION 

C 

c 

99  CONTINUE 
EnO 
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END 

DATE 

FILMED 


\ 


